---
title: "Beck Hall 2018 Experimental Analysis Code"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

library(lsmeans)
library(car)
library(ggplot2)
library(ggpubr)
library(doBy)
library(RColorBrewer)
```

##Data Set-Up

```{r}

##Raw Data Files for Plots and one-way ANOVAs comparing treatments to controls

DataTempRaw=read.csv("BeckHall_2018_Expts_RawData.csv")
Data=subset(DataTempRaw, Expt=="A"| Expt=="B") #Subset to include expts 1 and 2
Data3=subset(DataTempRaw, Expt=="C") #Subset to include expt 3
Data3$Conc=factor(Data3$Conc, levels=c("Control", "0.05 M", "0.5 M"))

#Log Response Ratio Data Files for one-way and two-way ANOVAs comparing treatment classes

DataTempLRR=read.csv("BeckHall_2018_Expts_LogResponseRatios.csv")
Data1=subset(DataTempLRR, Expt=="A"| Expt=="B") #Subset to include expts 1 and 2
Data2=subset(DataTempLRR, Expt=="C") #Subset to include expt 3
Data2$Conc=factor(Data2$Conc)

```

##One-way ANOVAs comparing all eight treatments to controls (Using Expts 1-2 Only)

###Chlorophyll a

```{r}

#Initial ANOVA
PFit=aov(Chla~Trt+Expt, data=Data)
summary(PFit)

#Testing for ANOVA Assumptions
par(mfrow=c(2,2))
plot(PFit)

#Dunnett
rg.Fit=ref.grid(PFit)
lsm.Fit= lsmeans(rg.Fit, "Trt")
comparison= contrast(lsm.Fit, "trt.vs.ctrl1")
comparison


```

###AFDM

```{r}

#Initial ANOVA
PFit=aov(AFDM~Trt+Expt, data=Data)
summary(PFit)

#Testing for ANOVA Assumptions
par(mfrow=c(2,2))
plot(PFit)

#Log ANOVA
PFit2=aov(log(AFDM)~Trt+Expt, data=Data)
summary(PFit2)

#Testing for ANOVA Assumptions
par(mfrow=c(2,2))
plot(PFit2)

#Dunnett
rg.Fit=ref.grid(PFit2)
lsm.Fit= lsmeans(rg.Fit, "Trt")
comparison= contrast(lsm.Fit, "trt.vs.ctrl1")
comparison

```

###Autotrophic Index

```{r}

#Initial ANOVA
PFit=aov(AI~Trt+Expt, data=Data)
summary(PFit)

#Testing for ANOVA Assumptions
par(mfrow=c(2,2))
plot(PFit)

PFit2=aov(log(AI)~Trt+Expt, data=Data)
summary(PFit)

#Testing for ANOVA Assumptions
par(mfrow=c(2,2))
plot(PFit2)

#Dunnett
rg.Fit=ref.grid(PFit2)
lsm.Fit= lsmeans(rg.Fit, "Trt")
comparison= contrast(lsm.Fit, "trt.vs.ctrl1")
comparison


```

###GPP

```{r}

#Initial ANOVA
PFit=aov(GPP~Trt+Expt, data=Data)
summary(PFit)

#Testing for ANOVA Assumptions
par(mfrow=c(2,2))
plot(PFit)

#Dunnett
rg.Fit=ref.grid(PFit)
lsm.Fit= lsmeans(rg.Fit, "Trt")
comparison= contrast(lsm.Fit, "trt.vs.ctrl1")
comparison


```

###GPP/Chla

```{r}
#Initial ANOVA
PFit=aov(GPP.Chla~Trt+Expt, data=Data)
summary(PFit)

#Testing for ANOVA Assumptions
par(mfrow=c(2,2))
plot(PFit)

#Dunnett
rg.Fit=ref.grid(PFit)
lsm.Fit= lsmeans(rg.Fit, "Trt")
comparison= contrast(lsm.Fit, "trt.vs.ctrl1")
comparison


```

###Respiration

```{r}

#Initial ANOVA
PFit=aov(Resp~Trt+Expt, data=Data)
summary(PFit)

#Testing for ANOVA Assumptions
par(mfrow=c(2,2))
plot(PFit)

#Dunnett
rg.Fit=ref.grid(PFit)
lsm.Fit= lsmeans(rg.Fit, "Trt")
comparison= contrast(lsm.Fit, "trt.vs.ctrl1")
comparison

```

##ANOVAs with Chlorophyll LRR as the Response Variable and Treatment Classes as Predictors

```{r}

##ANOVAs (Main effects, and Interactions with Concentration)

ChlaCat=aov(Chla~Cation+ Expt, data=Data1)
summary(ChlaCat)

#Testing for ANOVA Assumptions
par(mfrow=c(2,2))
plot(ChlaCat)

ChlaHyd=aov(Chla~Hydrogens+ Expt, data=Data1)
summary(ChlaHyd)

ChlaHeat=aov(Chla~Heat+ Expt, data=Data1)
summary(ChlaHeat)

options(contrasts=c("contr.sum", "contr.poly"))
ChlaCatConc=aov(Chla~Cation*Conc, data=Data2)
Anova(ChlaCatConc, type="III") 

options(contrasts=c("contr.sum", "contr.poly"))
ChlaHydConc=aov(Chla~Hydrogens*Conc, data=Data2)
Anova(ChlaHydConc, type="III") 

options(contrasts=c("contr.sum", "contr.poly"))
ChlaHeatConc=aov(Chla~Heat*Conc, data=Data2)
Anova(ChlaHeatConc, type="III") 

```

##ANOVAs with AFDM LRR as the Response Variable and Treatment Classes as Predictors

```{r}

##ANOVAs (Main effects, and Interactions with Concentration)

AFDMCat=aov(AFDM~Cation+ Expt, data=Data1)
summary(AFDMCat)

#Testing for ANOVA Assumptions
par(mfrow=c(2,2))
plot(AFDMCat)

AFDMHyd=aov(AFDM~Hydrogens+ Expt, data=Data1)
summary(AFDMHyd)

AFDMHeat=aov(AFDM~Heat+ Expt, data=Data1)
summary(AFDMHeat)

options(contrasts=c("contr.sum", "contr.poly"))
AFDMCatConc=aov(AFDM~Cation*Conc, data=Data2)
Anova(AFDMCatConc, type="III") 

options(contrasts=c("contr.sum", "contr.poly"))
AFDMHydConc=aov(AFDM~Hydrogens*Conc, data=Data2)
Anova(AFDMHydConc, type="III") 

lsmeans(AFDMHydConc, pairwise~Hydrogens)

options(contrasts=c("contr.sum", "contr.poly"))
AFDMHeatConc=aov(AFDM~Heat*Conc, data=Data2)
Anova(AFDMHeatConc, type="III") 

```

##ANOVAs with AI LRR as the Response Variable and Treatment Classes as Predictors

```{r}

##ANOVAs (Main effects, and Interactions with Concentration)

AICat=aov(AI~Cation+ Expt, data=Data1)
summary(AICat)

#Testing for ANOVA Assumptions
par(mfrow=c(2,2))
plot(AICat)

AIHyd=aov(AI~Hydrogens+ Expt, data=Data1)
summary(AIHyd)

AIHeat=aov(AI~Heat+ Expt, data=Data1)
summary(AIHeat)

options(contrasts=c("contr.sum", "contr.poly"))
AICatConc=aov(AI~Cation*Conc, data=Data2)
Anova(AICatConc, type="III") 

options(contrasts=c("contr.sum", "contr.poly"))
AIHydConc=aov(AI~Hydrogens*Conc, data=Data2)
Anova(AIHydConc, type="III") 

options(contrasts=c("contr.sum", "contr.poly"))
AIHeatConc=aov(AI~Heat*Conc, data=Data2)
Anova(AIHeatConc, type="III") 

```

##ANOVAs with GPP LRR as the Response Variable and Treatment Classes as Predictors

```{r}

##ANOVAs (Main effects, and Interactions with Concentration)

GPPCat=aov(GPP~Cation+ Expt, data=Data1)
summary(GPPCat)

#Testing for ANOVA Assumptions
par(mfrow=c(2,2))
plot(GPPCat)

GPPHyd=aov(GPP~Hydrogens+ Expt, data=Data1)
summary(GPPHyd)

lsmeans(GPPHyd, pairwise~Hydrogens)

GPPHeat=aov(GPP~Heat+ Expt, data=Data1)
summary(GPPHeat)

options(contrasts=c("contr.sum", "contr.poly"))
GPPCatConc=aov(GPP~Cation*Conc, data=Data2)
Anova(GPPCatConc, type="III") 

options(contrasts=c("contr.sum", "contr.poly"))
GPPHydConc=aov(GPP~Hydrogens*Conc, data=Data2)
Anova(GPPHydConc, type="III") 

options(contrasts=c("contr.sum", "contr.poly"))
GPPHeatConc=aov(GPP~Heat*Conc, data=Data2)
Anova(GPPHeatConc, type="III") 

```

##ANOVAs with Resp LRR as the Response Variable and Treatment Classes as Predictors

```{r}

##ANOVAs (Main effects, and Interactions with Concentration)

RespCat=aov(Resp~Cation+ Expt, data=Data1)
summary(RespCat)

#Testing for ANOVA Assumptions
par(mfrow=c(2,2))
plot(RespCat)

lsmeans(RespCat, pairwise~Cation)

RespHyd=aov(Resp~Hydrogens+ Expt, data=Data1)
summary(RespHyd)

RespHeat=aov(Resp~Heat+ Expt, data=Data1)
summary(RespHeat)

options(contrasts=c("contr.sum", "contr.poly"))
RespCatConc=aov(Resp~Cation*Conc, data=Data2)
Anova(RespCatConc, type="III") 

options(contrasts=c("contr.sum", "contr.poly"))
RespHydConc=aov(Resp~Hydrogens*Conc, data=Data2)
Anova(RespHydConc, type="III") 

options(contrasts=c("contr.sum", "contr.poly"))
RespHeatConc=aov(Resp~Heat*Conc, data=Data2)
Anova(RespHeatConc, type="III") 

```


##ANOVAs with GPP/Chla LRR as the Response Variable and Treatment Classes as Predictors

```{r}

##ANOVAs (Main effects, and Interactions with Concentration)

GPP.ChlaCat=aov(GPP.Chla~Cation+ Expt, data=Data1)
summary(GPP.ChlaCat)

#Testing for ANOVA Assumptions
par(mfrow=c(2,2))
plot(GPP.ChlaCat)

GPP.ChlaHyd=aov(GPP.Chla~Hydrogens+ Expt, data=Data1)
summary(GPP.ChlaHyd)

lsmeans(GPP.ChlaHyd, pairwise~Hydrogens)

GPP.ChlaHeat=aov(GPP.Chla~Heat+ Expt, data=Data1)
summary(GPP.ChlaHeat)

options(contrasts=c("contr.sum", "contr.poly"))
GPP.ChlaCatConc=aov(GPP.Chla~Cation*Conc, data=Data2)
Anova(GPP.ChlaCatConc, type="III") 

options(contrasts=c("contr.sum", "contr.poly"))
GPP.ChlaHydConc=aov(GPP.Chla~Hydrogens*Conc, data=Data2)
Anova(GPP.ChlaHydConc, type="III") 

options(contrasts=c("contr.sum", "contr.poly"))
GPP.ChlaHeatConc=aov(GPP.Chla~Heat*Conc, data=Data2)
Anova(GPP.ChlaHeatConc, type="III")

lsmeans(GPP.ChlaHeatConc, pairwise~Heat)

```

##Barplots

```{r}

theme_set(theme_bw(base_size = 14))

sumfun <- function(x,...){
c(m=mean(x,...), l=length(x), stdev=sd(x,...))
}
Expt2Sum=summaryBy(Chla+AFDM+AI+GPP+GPP.Chla+Resp~Trt, data=Data, FUN=sumfun, na.rm=TRUE)
Expt2Sum$ChlaSE=Expt2Sum$Chla.stdev/(sqrt(Expt2Sum$Chla.l))
Expt2Sum$AFDMSE=Expt2Sum$AFDM.stdev/(sqrt(Expt2Sum$AFDM.l))
Expt2Sum$AISE=Expt2Sum$AI.stdev/(sqrt(Expt2Sum$AI.l))
Expt2Sum$GPPSE=Expt2Sum$GPP.stdev/(sqrt(Expt2Sum$GPP.l))
Expt2Sum$GPP.ChlaSE=Expt2Sum$GPP.Chla.stdev/(sqrt(Expt2Sum$GPP.Chla.l))
Expt2Sum$RespSE=Expt2Sum$Resp.stdev/(sqrt(Expt2Sum$Resp.l))

sumfun <- function(x,...){
c(m=mean(x,...), l=length(x), stdev=sd(x,...))
}
Expt3Sum=summaryBy(Chla+AFDM+AI+GPP+GPP.Chla+Resp~Conc, data=Data3, FUN=sumfun, na.rm=TRUE)
Expt3Sum$ChlaSE=Expt3Sum$Chla.stdev/(sqrt(Expt3Sum$Chla.l))
Expt3Sum$AFDMSE=Expt3Sum$AFDM.stdev/(sqrt(Expt3Sum$AFDM.l))
Expt3Sum$AISE=Expt3Sum$AI.stdev/(sqrt(Expt3Sum$AI.l))
Expt3Sum$GPPSE=Expt3Sum$GPP.stdev/(sqrt(Expt3Sum$GPP.l))
Expt3Sum$GPP.ChlaSE=Expt3Sum$GPP.Chla.stdev/(sqrt(Expt3Sum$GPP.Chla.l))
Expt3Sum$RespSE=Expt3Sum$Resp.stdev/(sqrt(Expt3Sum$Resp.l))

#Chlorophyll a

Chla1=ggplot(Expt2Sum, aes(x=Trt, y=Chla.m, fill=Trt)) +
  labs(x="", y=expression(Chlorophyll~paste(italic(a))~mu*g~cm^{-2})) +
  geom_bar(position=position_dodge(.9), width=1, colour="black", stat="identity") +
  geom_errorbar(position=position_dodge(.9), width=.25, aes(ymin=Chla.m-ChlaSE, ymax=Chla.m+ChlaSE))+
  #scale_fill_grey(start = 0.25, end = .9)+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(), panel.grid.minor = element_blank(), panel.grid.major = element_blank(),
        axis.text.y=element_text(colour="black"))+
      scale_fill_brewer(palette="OrRd", labels=c("Control", bquote("K"[2]*"H"*"-HS"), bquote("K"[2]*"H"*"-HT"), bquote("K"*"H"[2]*"-HS"), bquote("K"*"H"[2]*"-HT"), bquote("Na"[2]*"H"*"-HS"), bquote("Na"[2]*"H"*"-HT"), bquote("Na"*"H"[2]*"-HS"), bquote("Na"*"H"[2]*"-HT")))

Chla2=ggplot(Expt3Sum, aes(x=Conc, y=Chla.m, fill=Conc)) +
  labs(x="", y="") +
  geom_bar(position=position_dodge(.9), width=1, colour="black", stat="identity") +
  geom_errorbar(position=position_dodge(.9), width=.25, aes(ymin=Chla.m-ChlaSE, ymax=Chla.m+ChlaSE))+
  #scale_fill_grey(start = .25, end = .9)+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(), panel.grid.minor = element_blank(), panel.grid.major = element_blank(), axis.text.y=element_text(colour="black"))+
      scale_fill_brewer(palette="OrRd")

figure <- ggarrange(Chla1, Chla2, labels = c("A", "B"), ncol = 2, nrow = 1)
#annotate_figure(figure, 
#                top = text_grob("Chlorophyll a by treatment", face = "bold",hjust=.6, size = 18))


ggsave("Chla_plot.tiff", width=18, height=10, units="cm", dpi=300)

#AFDM

AFDM1=ggplot(Expt2Sum, aes(x=Trt, y=AFDM.m, fill=Trt)) +
  labs(x="", y=expression(AFDM~mu*g~cm^{-2})) +
  geom_bar(position=position_dodge(.9), width=1, colour="black", stat="identity") +
  geom_errorbar(position=position_dodge(.9), width=.25, aes(ymin=AFDM.m-AFDMSE, ymax=AFDM.m+AFDMSE))+
  #scale_fill_grey(start = .25, end = .9)+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(), panel.grid.minor = element_blank(), panel.grid.major = element_blank(), axis.text.y=element_text(colour="black"))+
  ylim(0,1500)+
      scale_fill_brewer(palette="OrRd", labels=c("Control", bquote("K"[2]*"H"*"-HS"), bquote("K"[2]*"H"*"-HT"), bquote("K"*"H"[2]*"-HS"), bquote("K"*"H"[2]*"-HT"), bquote("Na"[2]*"H"*"-HS"), bquote("Na"[2]*"H"*"-HT"), bquote("Na"*"H"[2]*"-HS"), bquote("Na"*"H"[2]*"-HT")))

AFDM2=ggplot(Expt3Sum, aes(x=Conc, y=AFDM.m, fill=Conc)) +
  labs(x="", y="") +
  geom_bar(position=position_dodge(.9), width=1,colour="black", stat="identity") +
  geom_errorbar(position=position_dodge(.9), width=.25, aes(ymin=AFDM.m-AFDMSE, ymax=AFDM.m+AFDMSE))+
  #scale_fill_grey(start = .25, end = .9)+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(), panel.grid.minor = element_blank(), panel.grid.major = element_blank(), axis.text.y=element_text(colour="black"))+
  ylim(0,2500)+
      scale_fill_brewer(palette="OrRd")

figure <- ggarrange(AFDM1, AFDM2,labels = c("A", "B"),
                    ncol = 2, nrow = 1)
#annotate_figure(figure,
 #               top = text_grob("AFDM by treatment", face = "bold", hjust=.6,size = 18))


ggsave("AFDM_plot.tiff", width=18, height=10, units="cm", dpi=300)


#AI

AI1=ggplot(Expt2Sum, aes(x=Trt, y=AI.m, fill=Trt)) +
  ylab("Autotrophic Index")+
  geom_bar(position=position_dodge(.9), width=1,colour="black", stat="identity") +
  geom_errorbar(position=position_dodge(.9), width=.25, aes(ymin=AI.m-AISE, ymax=AI.m+AISE))+
 # scale_fill_grey(start = .25, end = .9)+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(), panel.grid.minor = element_blank(), panel.grid.major = element_blank(), axis.text.y=element_text(colour="black"))+
      scale_fill_brewer(palette="OrRd", labels=c("Control", bquote("K"[2]*"H"*"-HS"), bquote("K"[2]*"H"*"-HT"), bquote("K"*"H"[2]*"-HS"), bquote("K"*"H"[2]*"-HT"), bquote("Na"[2]*"H"*"-HS"), bquote("Na"[2]*"H"*"-HT"), bquote("Na"*"H"[2]*"-HS"), bquote("Na"*"H"[2]*"-HT")))
  
AI2=ggplot(Expt3Sum, aes(x=Conc, y=AI.m, fill=Conc)) +
  ylab("")+
  geom_bar(position=position_dodge(.9),width=1, colour="black", stat="identity") +
  geom_errorbar(position=position_dodge(.9), width=.25, aes(ymin=AI.m-AISE, ymax=AI.m+AISE))+
  #scale_fill_grey(start = .25, end = .9)+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(), panel.grid.minor = element_blank(), panel.grid.major = element_blank(), axis.text.y=element_text(colour="black"))+
  ylim(0,6000)+
      scale_fill_brewer(palette="OrRd")

figure <- ggarrange(AI1, AI2,labels = c("A", "B"),
                    ncol = 2, nrow = 1)
#annotate_figure(figure,
 #               top = text_grob("Autotrophic index by treatment", face = "bold", size = 18))


ggsave("AI_plot.tiff", width=18, height=10, units="cm", dpi=300)


#GPP

GPP1=ggplot(Expt2Sum, aes(x=Trt, y=GPP.m, fill=Trt)) +
  labs(x="", y=expression(GPP~mu*g~O[2]~cm^{-2}~hr^{-1})) +
  geom_bar(position=position_dodge(.9), width=1,colour="black", stat="identity") +
  geom_errorbar(position=position_dodge(.9), width=.25, aes(ymin=GPP.m-GPPSE, ymax=GPP.m+GPPSE))+
  #scale_fill_grey(start = .25, end = .9)+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(), panel.grid.minor = element_blank(), panel.grid.major = element_blank(), axis.text.y=element_text(colour="black"))+
  ylim(0,80)+
      scale_fill_brewer(palette="OrRd", labels=c("Control", bquote("K"[2]*"H"*"-HS"), bquote("K"[2]*"H"*"-HT"), bquote("K"*"H"[2]*"-HS"), bquote("K"*"H"[2]*"-HT"), bquote("Na"[2]*"H"*"-HS"), bquote("Na"[2]*"H"*"-HT"), bquote("Na"*"H"[2]*"-HS"), bquote("Na"*"H"[2]*"-HT")))

GPP2=ggplot(Expt3Sum, aes(x=Conc, y=GPP.m, fill=Conc)) +
  labs(x="", y="") +
  geom_bar(position=position_dodge(.9),width=1, colour="black", stat="identity") +
  geom_errorbar(position=position_dodge(.9), width=.25, aes(ymin=GPP.m-GPPSE, ymax=GPP.m+GPPSE))+
  #scale_fill_grey(start = .25, end = .9)+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(), panel.grid.minor = element_blank(), panel.grid.major = element_blank(), axis.text.y=element_text(colour="black"))+
      scale_fill_brewer(palette="OrRd")

figure <- ggarrange(GPP1, GPP2, labels = c("A", "B"),
                    ncol = 2, nrow = 1)
#annotate_figure(figure,
  #              top = text_grob("GPP by treatment", hjust=.65, face = "bold", size = 18))


ggsave("GPP_plot.tiff", width=18, height=10, units="cm", dpi=300)


#Resp

Resp1=ggplot(Expt2Sum, aes(x=Trt, y=Resp.m, fill=Trt)) +
  labs(x="", y=expression(Respiration~mu*g~O[2]~cm^{-2}~hr^{-1})) +
  geom_bar(position=position_dodge(.9), width=1,colour="black", stat="identity") +
  geom_errorbar(position=position_dodge(.9), width=.25, aes(ymin=Resp.m-RespSE, ymax=Resp.m+RespSE))+
  #scale_fill_grey(start = .25, end = .9)+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(), panel.grid.minor = element_blank(), panel.grid.major = element_blank(), axis.text.y=element_text(colour="black"))+
      scale_fill_brewer(palette="OrRd", labels=c("Control", bquote("K"[2]*"H"*"-HS"), bquote("K"[2]*"H"*"-HT"), bquote("K"*"H"[2]*"-HS"), bquote("K"*"H"[2]*"-HT"), bquote("Na"[2]*"H"*"-HS"), bquote("Na"[2]*"H"*"-HT"), bquote("Na"*"H"[2]*"-HS"), bquote("Na"*"H"[2]*"-HT")))

Resp2=ggplot(Expt3Sum, aes(x=Conc, y=Resp.m, fill=Conc)) +
  labs(x="", y="") +
  geom_bar(position=position_dodge(.9),width=1, colour="black", stat="identity") +
  geom_errorbar(position=position_dodge(.9), width=.25, aes(ymin=Resp.m-RespSE, ymax=Resp.m+RespSE))+
  #scale_fill_grey(start = .25, end = .9)+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(), panel.grid.minor = element_blank(), panel.grid.major = element_blank(), axis.text.y=element_text(colour="black"))+
      scale_fill_brewer(palette="OrRd")

figure <- ggarrange(Resp1, Resp2, labels = c("A", "B"),
                    ncol = 2, nrow = 1)
#annotate_figure(figure,
  #              top = text_grob("Respiration by treatment", hjust=.6, face = "bold", size = 18))


ggsave("Resp_plot.tiff", width=18, height=10, units="cm", dpi=300)


#GPP/Chla

GPPC1=ggplot(Expt2Sum, aes(x=Trt, y=GPP.Chla.m, fill=Trt)) +
  labs(x="", y=expression(GPP/ Chlorophyll~paste(italic(a)))) +
  geom_bar(position=position_dodge(.9), width=1,colour="black", stat="identity") +
  geom_errorbar(position=position_dodge(.9), width=.25, aes(ymin=GPP.Chla.m-GPP.ChlaSE, ymax=GPP.Chla.m+GPP.ChlaSE))+
  #scale_fill_grey(start = .25, end = .9)+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(), panel.grid.minor = element_blank(), panel.grid.major = element_blank(), axis.text.y=element_text(colour="black"))+
  ylim(0,300)+
      scale_fill_brewer(palette="OrRd", labels=c("Control", bquote("K"[2]*"H"*"-HS"), bquote("K"[2]*"H"*"-HT"), bquote("K"*"H"[2]*"-HS"), bquote("K"*"H"[2]*"-HT"), bquote("Na"[2]*"H"*"-HS"), bquote("Na"[2]*"H"*"-HT"), bquote("Na"*"H"[2]*"-HS"), bquote("Na"*"H"[2]*"-HT")))

GPPC2=ggplot(Expt3Sum, aes(x=Conc, y=GPP.Chla.m, fill=Conc)) +
  labs(x="", y="")+
  geom_bar(position=position_dodge(.9), width=1,colour="black", stat="identity") +
  geom_errorbar(position=position_dodge(.9), width=.25, aes(ymin=GPP.Chla.m-GPP.ChlaSE, ymax=GPP.Chla.m+GPP.ChlaSE))+
  #scale_fill_grey(start = .25, end = .9)+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(), panel.grid.minor = element_blank(), panel.grid.major = element_blank(), axis.text.y=element_text(colour="black"))+
      scale_fill_brewer(palette="OrRd")

figure <- ggarrange(GPPC1, GPPC2,labels = c("A", "B"),
                    ncol = 2, nrow = 1)
#annotate_figure(figure,
    #            top = text_grob("Chlorophyll-specific GPP by treatment", face = "bold", size = 18))

ggsave("GPPChla_plot.tiff", width=18, height=10, units="cm", dpi=300)

```

##Lab Experiments

###pH of Agar

```{r}

Data=read.csv("BeckHall_2018_Expts_AgarpH.csv")
options(contrasts=c("contr.sum", "contr.poly"))
Mod=aov(pH~Hydrogens*Heat, data=Data)
Anova(Mod, type="III")

#Testing for ANOVA Assumptions
par(mfrow=c(2,2))
plot(Mod)

lsmeans(Mod, pairwise~Hydrogens)

```

###pH of Distilled Water

```{r}

Data=read.csv("BeckHall_2018_Expts_DIpH.csv")
options(contrasts=c("contr.sum", "contr.poly"))
Mod=aov(pH_Change~Hydrogens*Heat, data=Data)
Anova(Mod, type="III")

#Testing for ANOVA Assumptions
#Result: Does Meet assumptions
par(mfrow=c(2,2))
plot(Mod)

lsmeans(Mod, pairwise~Hydrogens)
```

